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Abstract: We find that a surface soliton in nonlocal nonlinear media can 
be regarded as a half of a bulk soliton with an antisymmetric amplitude 
distribution. The analytical solutions for the surface solitons and breathers 
in strongly nonlocal media are obtained, and the critical power and breather 
period are gotten analytically and confirmed by numerical simulations. In 
addition, the oscillating propagation of nonlocal surface solitons launched 
away from the stationary position is considered as the interaction between 
the soliton and its out-of-phase image beam. Its trajectory and oscillating 
period obtained by our model are in good agreement with the numerical 
simulations. 
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1. Introduction 

The nonlocality of the nonlinear response exists in many real physical systems, such as pho- 
torefractive crystals lH], nematic liquid crystals ill] |3] SI |5l , lead glasses 161121, atomic vapors 
Is), Bose-Einstein condensates mHOl etc. In nonlocal nonlinear media, various soliton solu- 
tions have been predicted theoretically, such as vortex solitons lITTl [121 . multi-pole solitons 
|17]lT3l[l4l[l5|, Laguerre-Gaussian and Hermite-Gaussian solitons |[T6l[r7HT8l . Ince-Gaussian 
solitons |[T9l . and some of these have been observed experimentally |l4ll5]l6ll2ll2Ql- For nonlocal 
solitons, there are many interesting properties, for instance, the large phase shift 1211 . attraction 
between two bright out-of-phase solitons 1511221, attraction between two dark solitons Il23ll24l . 
self-induced fractional Fourier transform [ESl . etc. 

Recently, the nonlocal surface solitons have been investigated numerically and experimen- 
tally ll26l l27l l28l |29l [30l ISTl . The nonlocal surface solitons occurring at the interface between 
a dielectric medium and a nonlocal material exhibit unique properties. Nonlocal multipole sur- 
face solitons, vortices, and bound states of vortex solitons, incoherent surface solitons and ring 
surface waves have been also studied recently lIZTl l28l l29l l30ll . But the properties of nonlocal 
surface solitons mentioned above are all discussed by numerical simulations. In order to get a 
good understanding of the properties of nonlocal surface solitons, it is essential to present an 
analytic solution, even an approximate one. 

In this paper, based on the propagation equations governing the nonlocal surface waves and 
assuming that all energy of the surface soliton resides in nonlocal media, we find a surface 
soliton in nonlocal nonlinear media can be regarded as a half of a bulk soliton with an anti- 
symmetric amplitude distribution. The evolution regularities for the nonlocal surface waves are 
discussed both analytically and numerically. The analytical solutions for the surface solitons 
and breathers in strongly nonlocal media are obtained, and the critical power and breather pe- 
riod are gotten analytically and confirmed by numerical simulations. In addition, the oscillating 
propagation of nonlocal surface solitons launched away from the stationary position is consid- 
ered as the interaction between the soliton and its out-of-phase image beam. Its trajectory and 
oscillating period obtained by our model are in good agreement with the numerical simulations. 

2. Relation between surface solitons and bulk solitons 



Nonlocal medium Linear medium 

x<0 I x>0 



Nonlocal bulk medium 

x<0 , x>0 



Fig. 1. Sketches of a nonlocal surface soliton (a) and a nonlocal antisymmetric bulk soliton 
(b). Solid and dashed lines represent the distributions of the amplitude and the nonlinear 
refraction index, respectively. 



We consider a (iH-l)-D model of an optical surface wave with an envelope q propagating in z 
direction near an interface between a nonlocal nonlinear medium and a linear medium [see Fig. 
[TJa)]. The propagation of surface waves is governed by the dimensionless nonlocal nonlinear 
Schrodinger equation (NNLSE), i.e. 



(i) in nonlocal nonlinear media, x < 0, 



(ii) in linear media, x > Q, 



dq 1 d^q 

/^ + -^+A«, = 0, (1) 
A«-w,i^ = \q\\ (2) 

where jc,z denote, respectively, the normalized transversal and longitudinal coordinates. A« is 
the nonlinear perturbation of the refractive index, and denotes the normalized difference 
between the unperturbed refractive index of the nonlocal medium and the refractive index of 
the less optically dense linear medium. w,„ represents the characteristic length of the nonlocal 
material response, and a = w,n/wQ is the degree of nonlocality, where wq is the beam width. 
The optical intensity is normalized as / = \q'\^. In fact, Eqs. ([TJ and (|2| denote a NNLSE with 
an exponential-decay type nonlocal response. The exponential-decay type nonlocal response 
exists in many real physical systems, for instance, all diffusion-type nonlinearity 1321 [33l . 
orientational-type nonlinearity ||2]|22], and the general quadratic nonlinearity describing para- 
metric interaction Il35ll36l . Moreover when w,„ — > oo, Eqs. ([T]i and Q can be transformed to the 
forms describing the thermal nonlocal media (for example, lead glassesll6l l26l ). 

The boundary conditions for the surface soliton are q{x — oo) = 0, q{+Q) = q{—Q), 
An{x-)- -oo) = 0, and dAn/ 3x^=0 = ||26l |27l ED . References 1261 E?] |2ll indicate if the 
index difference is big enough, namely n^/ ^ 1, which is easily satisfied in the actual physi- 
cal system, the optical energy is almost totally confined in the nonlocal medium. Therefore one 
can approximately get the relation q{—0) = ^(+0) = and q{x > 0) = 0. Under this approxi- 
mation, the propagation of surface waves can be solved totally only based on Eqs. ([U and ^ 
with the boundary conditions {x < 0), 

^(x^ -oc) =^(0) =0, (4a) 
dAn 

^ =0. (4b) 

We find the solution of a surface soliton under the approximation q{Q) = is identical with 
the half part of an antisymmetric soliton in a bulk medium as shown in Fig.[Tlb). For the bulk 
soliton solutions with antisymmetric amplitude distribution which is also governed by Eqs. ([T]i 
and one can easily obtain the following relations. 

^(-0)=<7(+0)=0, (5a) 
q{x ±oo) = 0, (5b) 
dAn 



dx 



= 0. (5c) 

x=0 



Comparing Eq. with Eq. (jSj, it can be found that the conditions that the surface solitons sat- 
isfy are the same as those of the left half (x < 0) of the antisymmetric bulk solitons. Therefore, 
the surface soliton can be regarded as a half of the bulk soliton, and all the results of antisym- 
metric bulk soliton can be transfered to the surface soliton. On the other hand, the right half of 
the bulk soliton can be regarded as an image beam of the surface soliton Il37i . The interaction 
between the surface soliton and the interface can be regarded as the interaction between the 
soliton and its image beam in bulk medium. 



In nonlocal bulk media, the soliton solutions can have the antisymmetric amplitude distribu- 
tion lfT3lll4l[T5]| . Some authors have demonstrated that Hermite-Gaussian function can be ap- 
plied to describe this type soliton solution in nonlinear media with several different nonlocal 
responses, especially for the strongly nonlocal case llTSi [T6l [iTl . In addition, it has been dis- 
covered that in a nonlinear material with a finite-range nonlocal or a very long-range nonlocal 
response, the maximal number of peaks in stable multipole bulk solitons is four, and all higher- 
order soliton bound states are unstable lfT3]fT4l . Therefore, surface solitons only with less than 
three poles can be stable in the case of thermally nonlocal interface, which is firstly addressed 
by Kartashov et a/ ||27l . 

In the following, we give some comparisons between the bulk solitons and the surface soli- 
tons to illustrate our conclusions, and some analytical results are given for the nonlocal surface 
solitons. 



3. Surface solitons and breathers 



In strongly nonlocal nonlinear media, there exist Hermite-Gaussian solitons or bre athers ifTSi [T6l 
[TTI . Because the first-order Hermite-Gaussian beam is antisymmetric about the beam center, 
it can be used to describe the nonlocal surface soliton. In bulk nonlocal media, the first-order 
Hermite-Gaussian trial beams can be expressed as 



q{x,z) — a(z)xexp 



w2(z) 



(6) 



where a(z),w(z),c(z), represent the amplitude, beam width, phase-front curvature and 
phase of the beams, respectively. The input power is fo = / Iq'P'ij^. 

In the following, Eqs. ^ and ^ are used for the equivalent bulk case (—00 < x < 0°). There- 
fore Eqs.([T]i and Q can be restated as an Euler-Lagrange equation corresponding to a varia- 
tional problem | 



where 



[L]dz = 0, 



\l ^ dq dq* dq dq* 
ax ax az az 



The Lagrangian density L is given by 
dq dq 



^ dz ^ <3z 



dq 
dx 



^-\q\^J^y{x~x')\q{x')\'dx', 



(7) 
(8) 

(9) 



where = (l/2w,„)exp(— |jic|/w„,). Because the complexity of integration produced by 

the expression of [L] is too difficult to obtain. But for the strongly nonlocal case, i.e. a = 

Wm/wo ^ 1, Rix) can be expanded as 



(10) 



2wi, 



2w2 



Combining Eq. (|6]l and Eq. ([SJ, [L] can be obtained 



[L] = --W-fl2H;(l+cV) -— W-flV ( 3w2— +4— 



.2., ,3 



512w 



2 ' ' 16 V 2 



dz dz 



(11) 



Following the common process of the variational approach, one can get four differential equa- 
tions, 

^ da ^ dw ^ , 

2w h3a — = 0, (12a) 

dz dz 

da ( dw \ 
2w— +fl 5— -4cvv =0, (12b) 
dz \ dz J 

2 5 

a w , 2 ~, /— -, 2\ 

— — (4;rw„, — l\/nwmW + inw ) 

- 4B.V27: (1 + c^w'^) -SVTjiw^ ( 3w^—+4-—] =0, (12c) 

\ dz dz J 

480 V2w,^c V'* + a^w^ {~24-y^wl + 49w,„w - 

+ 48V2w?„ f 4^2— +5w'*— + 2 ) =0. (12d) 
\ dz dz J 

Based on the above four equations, we can investigate the propagation of the trial beams. In the 
following, we divide the problem into two cases to discuss, i.e. solitons and breathers. 

3.1. surface solitons 

For the case of solitons, in Eq. (|6]l, a and w reduce to constants, and c = 0, 0(z) = j8z, where j3 
is propagation constant. The expression of beams can be rewritten in a simple form 

q{x,z) = aoxexp ( e'l^^. (13) 

Based on Eq. ( fT2] ). we can get the input power, namely the critical power of bulk solitons 

= , (14) 

wl(^w,n-6^/7iwo) ' 



and the propagation constant 

^_ 105w,„wo-6^/7i{8wl+9wl) ^^^^ 
2wl{6^wo - 7w,„) 

According to our conclusion in Sec. |2] the surface soliton can be expressed as 

( ( ,2 



1< 



, flo-xexp — e'*^' for x < 0, 
q{x,z)^{ V "'5/ - (16) 

for x>0. 



Apparently, the critical power of surface soltions is a half of bulk solitons, namely = Pc/2. 
When the optical beam propagates as a surface soliton, we define the position where the beam 
has its maximal intensity as the soliton position, i.e. the distance between the maximal intensity 
and the interface. According to Eq. ([16), the soliton position can be obtained as 

^/2^vo 

Xs = z — • (17) 



Based on Eqs. ([T|i-©i some numerical simulations of the propagations of surface and bulk 
solitons are carried out by using Eqs. (ITji and (IT6i as incident profiles. Figure|2]shows the ana- 
lytical and simulated results are in good agreement with each other in strongly nonlocal media. 
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Fig. 2. Comparison between the fundamental surface solitons and the first-order Hermite- 
Gaussian bulk solitons with the degree of nonlocality a = 10. (a) and (c) are simulated 
intensity distribution for surface solitons during propagation with wq = 1.0 and 2.0, re- 
spectively; (b) and (d) are simulated intensity distribution for bulk solitons with wq= 1.0 
and 2.0, respectively, (e)-(h) are, respectively, the transversal intensity distributions (solid 
line) and the refractive index distributions (dashed line) corresponding to row 1 . 
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Fig. 3. Simulated propagation of surface solitons in the presence of 3% white noise at 
a fixed nonlocal degree a = 10. (a) and (b) give the propagation of surface solitons with 
different soliton widths wq = 1.0 and 2.0, respectively, (c) and (d) give the incident intensity 
distribution (left) and the output intensity distribution at z = 2000 (right) corresponding to 
(a) and (b), respectively. 



In fig. 121 1 ~\q\^ is the optical intensity. To confirm the stabiHty, we simulate the propagation of 
surface solitons in the presence of 3% white noise (see fig.|3). Figure |3] shows the propagation 
is stable. 

In general nonlocal media, the profiles of surface soliton are obtained by numerical iterative 
method based on Eqs. ([Tli-(l3]i, and the relations between critical powers and beam width are 
shown in Fig. |4l For Figs. Ua) and|4jb) which belong to the strongly nonlocal case, the analyt- 
ical results are in excellent agreement with the numerical results. As the degree of nonlocality 
decreases to a = 6 in Fig. Etc), the analytical results are also accordant with the numerical 
results approximately. For Fig. SJd), the analytical results begin to deviate from the exact nu- 
merical ones with the beam width increasing. Figure 2] confirms that the analytical solutions 
of surface solitons are very valid for the strongly nonlocal case. The relations between criti- 
cal powers and propagation constants are shown in Fig. |5] It shows the analytical results from 
Eqs. ( fT4l i and ( fTSl l are in good agreement with the numerical results except for small j3 which 
corresponds to a weak nonlocality. 




Fig. 4. Variations of the surface soliton power with the soliton width under different condi- 
tions, (a) and (c) are the results with a fixed degree of nonlocality a = 30,6, respectively, 
(b) and (d) are the results with a fixed width of the nonlocal response w,„ = 100, 10, respec- 
tively. Solid lines represent the analytical results based on Eq. l ll4t : the triangles represent 
the numerical results. 
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Fig. 5. Variations of the surface soliton power with the propagation constant for (a) w„, = 30 
and (b) Wm = 10. Solid lines represent the analytical results based on Eqs. il4\ and il5\ \ 
the triangles represent the numerical results. 

For weakly nonlocal case, the analytical solution can not be obtained because the analytical 
solution for multipole bulk solitons can not be found in weakly nonlocal medium. However, in 
nonlocal bulk media, the diploe solitons are always existent and stable in the entire domain of 



their existence llT3l . According to the theory in Sec. |2] the fundamental surface soliton should 
be existent for the weakly nonlocal case. Following the method in Ref. lfT3]| . we seek the dipole 
bulk soliton combining two out-of -phase fundamental solitons under weakly nonlocality, as a 
result, the surface soliton can be also found at the same time, as shown in Fig.|6] 
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Fig. 6. Comparison between the surface solitons and the bulk solitons for the weakly nonlo- 
cal case (a = 0.2). (a) and (d) are the simulated intensity distributions for the surface soliton 
and the bulk soliton during propagation with = 1.0, vv„, = 0.5, respectively; (b) and (e) 
are the transversal intensity distributions (solid line) and the refractive index distributions 
(dashed line) corresponding to (a) and (d), respectively, (c) and (f) are the transversal ampli- 
tude distributions (dashed-dotted line) and the refractive index distributions (dashed line) 
corresponding to (a) and (d), respectively. 



3.2. surface breathers 

Because the breathers exist in nonlocal bulk media, if the input power is not equal to the soliton 
power, we can find the surface breathers. According to Eq. (IT2] ). we can get 



2V2 laW ^fna^w^ 



96 Vlwl 



I6wi 



(18) 



Above equation is equivalent to Newton's second law in classical mechanics for the motion of 
an one-dimensional particle with the equivalent mass 1 , thus one can get the equivalent potential 



y(w) 



V2 

2 



'TtW 



7a w' 



80wl 384-wl 



-co, 



(19) 



where cq is a constant. Through expanding the equivalent potential V{w) into the second order 
at the balance position, the motion of the particle is approximated to a harmonic oscillation, 
and we can obtain the approximate expression of the beam width 



V2 
2 



g — Vb + {V2wQ —g + Vb) cos 




12 



1 



(20) 



where 



lw„ 



b 
t 



49w2 



49wow,^, 

« = , 

1927r ' 

_ w§Pc(6y7rwo-7w,„) 

^ ~ 18^/7rPo 

According to our theory, the oscillation period of surface breathers is the same as that of bulk 
breathers. Then, we can obtain the oscillation period of the surface breather, 



Az = 2% 



12 



wlPs{lw,n-6^WQ) 



(21) 



where P^q is the input power of the surface breathers. At the same time, the beam trajectory can 
be obtained and given by xi, — —^/2w{z) /2. Figure|2]gives the comparison between the surface 
breathers and the bulk breathers for a = 10. The analytical trajectories are in good agreement 
with the simulated results. 




Fig. 7. Comparison between the surface breathers and the bulk breathers with a = 10. 
(a) and (c) are, respectively, the simulated propagations of the surface breather and the 
bulk breather with input power being twice soliton power, (b) and (d) give the analytical 
trajectories of the maximal intensity corresponding to (a) and (c), respectively. Row 2 is 
the same as row 1 except that the input power is a half of soliton power. 



4. Surface solitons launched away from the stationary position 

One can learn from Ref. ||26l that if the beam launched away from the soliton position, the 
beam will propagate oscillating about the stationary position. The oscillation induced by the 
interaction between the soliton and the interface is periodic, and the beam never converges 
to a straight line trajectory. According to our theory and introducing an out-of-phase image 



beam ll37l . the interaction between soliton and interface can be regarded as the interaction 
between the soliton and its image beam in nonlocal bulk media. Then the oscillating trajectory 
of surface solitons can be obtained from the trajectories of two out-of -phase solitons interacting 
in nonlocal bulk media. 

For the strongly nonlocal case, the fields of two out-of -phase Gaussian solitons propagating 
in bulk media can be expressed as 



q±(x,z) = ifloexp 



± iu{z)x 



(22) 



where 9{z) is the phase of whole beams, ±Xc{z) are the mass centers of Gaussian solitons and 
Xc(0) = XcQ- u{z) is the tilting wavefront induced by the attraction between two solitons, and we 
consider the normally incident case, i.e. m(0) = 0. The amplitude Oq = 32V2w^/wQ[4wm(wQ + 
2Wm) ~ V^wo{wq +4w,^j)] is normalized to guarantee the incident power of each soliton is 
equal to the critical power, so we can assume that the beam width does not change during 
propagation. 

If the two solitons are separated over 3wo each other, their fields almost do not overlap. 
However because of strong nonlocality, there still exists attractive force between them. The 
attractive force exerted on one soliton is produced by the other soliton, specifically, by the 
nonlinear refractive index change induced by the other soliton. Therefore, the trajectory of one 
soHton (i.e. the left one q+) can be determined by the light ray equation. 



dx 



= An', 



(23) 



where the refractive index change A«_ is only induced by the right soliton q^, as 



Rix-x')\q^{x',z)\^dx'. 



(24) 



Using the nonlinear response function R{x) = (1 /2vv„i)exp(— |jc|/wm) and Eqs. (I22li-(l24li. one 
can obtain 



d Xc 
'd? 



4V2w2 



-exp 



exp — erfc 



iWinX^ 



2V2i 



erf 



2^2' 



WQW„ 



1 



(25) 



where erf( ) and erfc( ) denote the error and complementary error functions, respectively. Tak- 
ing the conditions vv„, ^ wq^Xcq, we can get a simply approximate expression 



d^Xc 
'd^ 



2V2i 



(26) 



An' shown in Fig.[8le) represents the attractive force between two solitons. Due to the compli- 
cated form of Isii' , Eqs. (|25l l and (|26] l are solved numerically using Runge-Kutta methods. The 
soliton trajectories and the oscillation period are shown in Figs.[8jc) and[8jd). The trajectory of 
the surface soliton corresponds to the left half of trajectories of two solitons in bulk media, i.e. 

^sc — \-^c\- 

The numerical simulations based on Eqs. ([T]i - (|3]l for surface soliton launched away from 
stationary position and for two bulk solitons are also shown in Figs.[8ta) and[8tb). The trajecto- 
ries and the oscillation periods gotten from numerical simulations are compared with that from 



Fig. 8. (a) Simulated propagations of a Gaussian surface wave launched 2.5w() away from 
the surface, (b) Simulated propagations of two out-of-phase Gaussian solitons with a sep- 
arated distance S.Owq. (c) The trajectories of surface wave and two Gaussian solitons. (d) 
The oscillation period of surface wave, (e) Distribution of An'. The red solid line and the 
blue dashed-dotted line are the results based on Eqs. ( 125) and i26l . respectively. The green 
dashed line and green triangles are the simulated results directly based on Eqs. ([T}-(|3}. All 
cases are with the parameters Wm = 50.0, wq = 1.0. 

Eqs. ( IZSl l and ( |26] l in Figs.|8jc) and[8ld). These results are approximately coincident with each 
other. The main reason of disagreement is that Eq. (l2Jt is valid when two beam fields do not 
overlap. When Xco is small, the fields of two beams overlap and we can not distinguish A«+ and 
A«_ . Then Eq. (|23T l is invalid, and the simulated results disagree with that obtained by solving 
Eqs.(|25]l and ( [26] l. Whenx^o is large and the nonlocality is strong enough, the simulated results 
are in good agreement with the results by numerically solving Eqs. (|25] l and ( |26] l. 

5. Conclusion 

In conclusion, we have studied the nonlocal surface waves numerically and analytically. We find 
a surface soliton in nonlocal nonlinear media can be regarded as a half of a bulk soliton with an 
antisymmetric amplitude distribution. By applying the variational method and taking the first- 
order Hermite-Gaussian beam as an example, the analytical solutions for the surface solitons 
and breathers in strongly nonlocal media are obtained, and the critical power and breather pe- 
riod are gotten analytically and confirmed by numerical simulations. In addition, the oscillating 
propagation of nonlocal surface soliton launched away from the stationary position is consid- 
ered as the interaction between the soliton and its out-of-phase image beam. We have discussed 
the oscillation period and the beam trajectory. Its trajectory and oscillating period obtained by 
our model are in good agreement with the numerical simulations. 
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